# Making univariate toy data
library(dplyr,
quietly = TRUE)
mz_data <- MASS::mvrnorm(
n = 100,
mu = c(0, 0),
Sigma = matrix(c(1.0, 0.8,
0.8, 1.0),
nrow = 2)) %>%
as.data.frame() %>%
setNames(c('X1', 'X2'))
dz_data <- MASS::mvrnorm(
n = 100,
mu = c(0, 0),
Sigma = matrix(c(1.0, 0.5,
0.5, 1.0),
nrow = 2)) %>%
as.data.frame() %>%
setNames(c('X1', 'X2'))
data <- rbind(data.frame(zyg = 'MZ', mz_data),
data.frame(zyg = 'DZ', dz_data)) %>%
mutate(
zyg = factor(zyg, levels = c('MZ', 'DZ'))
)
head(data)
## zyg X1 X2
## 1 MZ 1.97640719 0.8614917
## 2 MZ 0.19813515 0.5722313
## 3 MZ -0.21476752 -0.1499316
## 4 MZ -0.05976396 -1.6416417
## 5 MZ -0.10528725 -0.8397834
## 6 MZ 1.39832330 1.1882021# Installing the packages
# install.packages('devtools')
# library(devtools)
# install_github('ivanvoronin/mlth.data.frame')
# install_github('ivanvoronin/TwinAnalysis')
# Loading the packages
library(OpenMx)
library(TwinAnalysis)
# Descriptive statistics
get_univ_descriptives(data, zyg = 'zyg', vars = 'X')
## $X
## Twin1------------------ Twin2------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 0.2892160 0.996696 100 0.1193898 0.8943036 0.7601622 0.6626901 0.8323086
## DZ 100 -0.1712873 1.067547 100 -0.1975880 0.9445352 0.4553250 0.2843362 0.5982401# Univariate ACE model
model <- univariate_ace(data, zyg = 'zyg', ph = 'X')
model <- mxRun(model)
summary(model)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -0.008753464 0.06321116
## 2 ACE.a[1,1] a X X 0.785090384 0.10184093
## 3 ACE.c[1,1] c X X 0.381478731 0.19748414
## 4 ACE.e[1,1] e X X 0.482854666 0.03427370
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 4 396
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1021.598
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/400
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 229.5979 1029.598
## BIC: -1076.5358 1042.791
## | Sample-Size Adjusted
## AIC: 1029.803
## BIC: 1030.119
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:43:25
## Wall clock time: 0.07648754 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Constrained models
# This is a list of MxModels
nested_models <- def_nested_models(
model,
AE = mxConstraint(c[1, 1] == 0),
CE = mxConstraint(a[1, 1] == 0),
E = list(mxConstraint(c[1, 1] == 0),
mxConstraint(a[1, 1] == 0)),
run = TRUE
)
fit_stats(model, nested_models = nested_models)
## Warning in computeFitStatistics(likelihood, DoF, chi,
## chiDoF, retval[["numObs"]], : Your model may be mis-
## specified (and fit worse than an independence model),
## or you may be using the wrong independence model, see ?
## mxRefModels
## model base ep minus2LL df AIC BIC
## 1 ACE Saturated 4 1021.598 396 229.5979 -1076.5358
## 2 AE ACE 4 1022.470 397 228.4697 -1080.9623
## 3 CE ACE 4 1039.653 397 245.6534 -1063.7786
## 4 E ACE 4 1130.385 398 334.3846 -978.3457
## CFI TLI RMSEA diffLL diffdf
## 1 0.8577940 0.9525980 0.1128539 21.2832071 6
## 2 0.8589871 0.9597106 0.1040432 0.8717727 1
## 3 0.6990974 0.9140278 0.1519838 18.0554872 1
## 4 -0.1358265 0.7160434 0.2762131 108.7866932 2
## p
## 1 1.631543e-03
## 2 3.504650e-01
## 3 2.145594e-05
## 4 2.383799e-24# Making bivariate toy data
mz_data2 <- MASS::mvrnorm(
n = 100,
mu = c(0, 0, 0, 0),
Sigma = matrix(
c(1.0, 0.7, 0.6, 0.5,
0.7, 1.0, 0.5, 0.8,
0.6, 0.5, 1.0, 0.7,
0.5, 0.8, 0.7, 1.0),
nrow = 4)) %>%
as.data.frame() %>%
setNames(c('X1', 'Y1', 'X2', 'Y2'))
dz_data2 <- MASS::mvrnorm(
n = 100,
mu = c(0, 0, 0, 0),
Sigma = matrix(
c(1.0, 0.7, 0.4, 0.2,
0.7, 1.0, 0.2, 0.4,
0.4, 0.2, 1.0, 0.7,
0.2, 0.4, 0.7, 1.0),
nrow = 4)) %>%
as.data.frame() %>%
setNames(c('X1', 'Y1', 'X2', 'Y2'))
data2 <- rbind(data.frame(zyg = 'MZ', mz_data2),
data.frame(zyg = 'DZ', dz_data2)) %>%
mutate(
zyg = factor(zyg, levels = c('MZ', 'DZ'))
)
head(data2)
## zyg X1 Y1 X2 Y2
## 1 MZ -0.5435463 -0.41234556 -0.86241804 -0.29155328
## 2 MZ 0.3948285 -0.84486187 -0.18510648 -1.11259359
## 3 MZ 0.3981069 -0.08203682 -0.05917079 -0.18544979
## 4 MZ -0.6577894 -0.73082395 -0.02673546 0.03067235
## 5 MZ 0.4123762 0.71153558 0.41007149 0.53892970
## 6 MZ -0.1671188 -1.43100131 -0.50763096 -1.31595716# Descriptive statistics
get_univ_descriptives(data2,
zyg = 'zyg',
vars = c('X', 'Y'))
## $X
## Twin1------------------- Twin2-------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 0.1602093 0.9891965 100 0.09714586 0.9135136 0.5967201 0.4534493 0.7099298
## DZ 100 -0.1450704 0.9265116 100 -0.20180809 1.0229349 0.4750329 0.3072861 0.6141479
##
## $Y
## Twin1------------------ Twin2-------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 0.1208040 1.101084 100 0.06408469 1.0693377 0.8602316 0.7987792 0.9039190
## DZ 100 -0.1001922 0.983983 100 -0.24580522 0.9595227 0.4185444 0.2420229 0.5682471# Cross-twin cross-trait correlations
get_cross_trait_cors(data2,
zyg = 'zyg',
vars = c('X', 'Y'))
## zyg: MZ
## X1 Y1 X2 Y2
## X1 1.0000000 0.6625874 0.5967201 0.4859931
## Y1 0.6625874 1.0000000 0.5098315 0.8602316
## X2 0.5967201 0.5098315 1.0000000 0.6378926
## Y2 0.4859931 0.8602316 0.6378926 1.0000000
## --------------------------------------------
## zyg: DZ
## X1 Y1 X2 Y2
## X1 1.0000000 0.6686251 0.4750329 0.2190231
## Y1 0.6686251 1.0000000 0.2353515 0.4185444
## X2 0.4750329 0.2353515 1.0000000 0.6493918
## Y2 0.2190231 0.4185444 0.6493918 1.0000000# Bivariate ACE model
model2 <- multivariate_ace(
data = data2,
zyg = 'zyg',
ph = c('X', 'Y')
)
model2 <- mxRun(model2)
summary(model2)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -3.293302e-02 0.06011835
## 2 ACE.mean[1,2] mean 1 2 -5.684037e-02 0.06549453
## 3 ACE.a[1,1] a X X 2.041926e-07 0.24116291
## 4 ACE.a[2,1] a Y X 5.774286e-01 0.10958100
## 5 ACE.a[2,2] a Y Y 9.063612e-01 0.09302743
## 6 ACE.c[1,1] c X X 4.986504e-01 0.24273142
## 7 ACE.c[2,1] c Y X -9.775818e-02 0.63518951
## 8 ACE.c[2,2] c Y Y 2.565577e-01 0.32877565
## 9 ACE.e[1,1] e X X 4.475017e-01 0.02362679
## 10 ACE.e[2,1] e Y X 3.870517e-01 0.04866025
## 11 ACE.e[2,2] e Y Y 4.042446e-01 0.02815001
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 11 789
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1793.14
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/800
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 215.1395 1815.140
## BIC: -2387.2329 1851.421
## | Sample-Size Adjusted
## AIC: 1816.544
## BIC: 1816.572
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:43:25
## Wall clock time: 0.07370257 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Output tables from the bivariate ACE
get_output_tables(model2)
## $Variance_components
## A(%) C(%) E(%)
## X 0.3540660 0.27419454 0.3717394
## Y 0.7818313 0.06264417 0.1555245
##
## $Covariation_components
## Total A C E
## X <-> X 0.9416995 0.3334238 0.25820885 0.3500668
## X <-> Y 0.6547418 0.5233589 -0.02508062 0.1564636
## Y <-> Y 1.0507262 0.8214906 0.06582187 0.1634137
## A(%) C(%) E(%)
## X <-> X 0.3540660 0.27419454 0.3717394
## X <-> Y 0.7424551 0.03558023 0.2219647
## Y <-> Y 0.7818313 0.06264417 0.1555245
##
## $All_correlations
## r_A r_C r_E Total
## X <-> X 1 1.0000000 1.0000000 1.0000000
## X <-> Y 1 -0.1923834 0.6541744 0.6582171
## Y <-> Y 1 1.0000000 1.0000000 1.0000000A model (mxModel):
# Twin models in OpenMx
library(OpenMx)
# Univariate ACE model
model <- mxModel(
name = 'ACE',
# Means
mxMatrix(type = 'Full',
nrow = 1, ncol = 1,
free = TRUE,
name = 'mean'),
mxAlgebra(cbind(mean, mean),
name = 'expMeans'),
# Variance components
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'a'),
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'c'),
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'e'),
mxAlgebra(t(a) %*% a, name = 'A'),
mxAlgebra(t(c) %*% c, name = 'C'),
mxAlgebra(t(e) %*% e, name = 'E'),
# Expected covariance
mxAlgebra(rbind(cbind(A + C + E, A + C ),
cbind(A + C, A + C + E)),
name = 'expCovMZ'),
mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
cbind(.5%x%A + C, A + C + E)),
name = 'expCovDZ'),
# MZ submodel
mxModel(name = 'MZ',
mxData(observed = mz_data,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovMZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'X2')),
mxFitFunctionML()),
# DZ submodel
mxModel(name = 'DZ',
mxData(observed = dz_data,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovDZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'X2')),
mxFitFunctionML()),
mxFitFunctionMultigroup(c('MZ','DZ')),
# Output tables
mxAlgebra(A + C + E, name = 'V'),
mxAlgebra(cbind(V, A, C, E),
dimnames = list('X', c('Total', 'A', 'C', 'E')),
name = 'Raw_variance'),
mxAlgebra(cbind(A / V, C / V, E / V),
dimnames = list('X', c('A(%)', 'C(%)', 'E(%)')),
name = 'Variance_components')
)
model <- mxRun(model)
summary(model)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -0.008753687 0.06321119
## 2 ACE.a[1,1] a 1 1 0.785090556 0.10183866
## 3 ACE.c[1,1] c 1 1 0.381478510 0.19747930
## 4 ACE.e[1,1] e 1 1 0.482854720 0.03427365
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 4 396
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1021.598
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/400
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 229.5979 1029.598
## BIC: -1076.5358 1042.791
## | Sample-Size Adjusted
## AIC: 1029.803
## BIC: 1030.119
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:43:25
## Wall clock time: 0.04036498 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Bivariate ACE model
model2 <- mxModel(
name = 'ACE',
# Means
mxMatrix(type = 'Full',
nrow = 1, ncol = 2,
free = TRUE,
name = 'mean'),
mxAlgebra(cbind(mean, mean),
name = 'expMeans'),
# Variance components
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'a'),
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'c'),
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'e'),
mxAlgebra(t(a) %*% a, name = 'A'),
mxAlgebra(t(c) %*% c, name = 'C'),
mxAlgebra(t(e) %*% e, name = 'E'),
# Covariance
mxAlgebra(rbind(cbind(A + C + E, A + C ),
cbind(A + C, A + C + E)),
name = 'expCovMZ'),
mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
cbind(.5%x%A + C, A + C + E)),
name = 'expCovDZ'),
# MZ submodel
mxModel(name = 'MZ',
mxData(observed = mz_data2,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovMZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'Y1', 'X2', 'Y2')),
mxFitFunctionML()),
# DZ submodel
mxModel(name = 'DZ',
mxData(observed = dz_data2,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovDZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'Y1', 'X2', 'Y2')),
mxFitFunctionML()),
mxFitFunctionMultigroup(c('MZ','DZ')),
# Output tables
mxAlgebra(A + C + E, name = 'V'),
mxAlgebra(abs(A) + abs(C) + abs(E), name = 'Vabs'),
mxAlgebra(abs(A) / Vabs, name = 'Aperc'),
mxAlgebra(abs(C) / Vabs, name = 'Cperc'),
mxAlgebra(abs(E) / Vabs, name = 'Eperc'),
mxAlgebra(cbind(diag2vec(Aperc),
diag2vec(Cperc),
diag2vec(Eperc)),
dimnames = list(c('X', 'Y'), c('A(%)', 'C(%)', 'E(%)')),
name = 'Variance_components'),
mxAlgebra(cbind(vech(V),
vech(A), vech(C), vech(E),
vech(Aperc), vech(Cperc), vech(Eperc)),
dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
c('Total', 'A', 'C', 'E', 'A(%)', 'C(%)', 'E(%)')),
name = 'Covariation_components'),
mxMatrix(type = "Iden",
nrow = 2, ncol = 2,
name = "I"),
mxAlgebra(sqrt(I * V), name = "SD_total"),
mxAlgebra(sqrt(I * A), name = 'SD_A'),
mxAlgebra(sqrt(I * C), name = 'SD_C'),
mxAlgebra(sqrt(I * E), name = 'SD_E'),
mxAlgebra(solve(SD_total) %&% V, name='r_total'),
mxAlgebra(solve(SD_A) %&% A, name = 'r_A'),
mxAlgebra(solve(SD_C) %&% C, name = 'r_C'),
mxAlgebra(solve(SD_E) %&% E, name = 'r_E'),
mxAlgebra(cbind(vech(r_A), vech(r_C), vech(r_E),
vech(r_total)),
dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
c('r_A', 'r_C', 'r_E', 'Total')),
name='All_correlations')
)
model2 <- mxRun(model2)
summary(model2)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -3.293282e-02 0.06011839
## 2 ACE.mean[1,2] mean 1 2 -5.684027e-02 0.06549455
## 3 ACE.a[1,1] a 1 1 -1.357998e-07 0.24112224
## 4 ACE.a[2,1] a 2 1 5.774276e-01 0.10955114
## 5 ACE.a[2,2] a 2 2 9.063608e-01 0.09299972
## 6 ACE.c[1,1] c 1 1 4.986521e-01 0.24264231
## 7 ACE.c[2,1] c 2 1 -9.775422e-02 0.63495296
## 8 ACE.c[2,2] c 2 2 2.565589e-01 0.32867213
## 9 ACE.e[1,1] e 1 1 4.475017e-01 0.02362661
## 10 ACE.e[2,1] e 2 1 3.870517e-01 0.04865963
## 11 ACE.e[2,2] e 2 2 4.042445e-01 0.02814985
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 11 789
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1793.14
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/800
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 215.1395 1815.140
## BIC: -2387.2329 1851.421
## | Sample-Size Adjusted
## AIC: 1816.544
## BIC: 1816.572
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:43:26
## Wall clock time: 0.06721616 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)mxEval(Covariation_components, model2)
## Total A C E
## X <-> X 0.9416993 0.3334227 0.25820980 0.3500668
## X <-> Y 0.6547416 0.5233578 -0.02507971 0.1564635
## Y <-> Y 1.0507260 0.8214899 0.06582245 0.1634136
## A(%) C(%) E(%)
## X <-> X 0.3540649 0.27419561 0.3717395
## X <-> Y 0.7424557 0.03557906 0.2219653
## Y <-> Y 0.7818308 0.06264473 0.1555245library(mlth.data.frame)
# mlth.data.frame - multiheaded data.frame
model2 <- multivariate_ace(
data = data2,
zyg = 'zyg',
ph = c('X', 'Y')
)
model2 <- def_ci(model2, ci = 'Variance_components')
model2 <- mxRun(model2, intervals = TRUE)
get_output_tables(model2)
## $Variance_components
## A(%)------------------------- C(%)----------------------------- E(%)-------------------------
## est lCI uCI est lCI uCI est lCI uCI
## X 0.3540660 0.1388926 0.5823833 0.27419454 9.228158e-02 0.4861962 0.3717394 0.2892442 0.4715775
## Y 0.7818313 0.5048375 0.8786609 0.06264417 6.563101e-31 0.3347188 0.1555245 0.1145364 0.2128753
##
## $Covariation_components
## Total A C E
## X <-> X 0.9416995 0.3334238 0.25820885 0.3500668
## X <-> Y 0.6547418 0.5233589 -0.02508062 0.1564636
## Y <-> Y 1.0507262 0.8214906 0.06582187 0.1634137
## A(%) C(%) E(%)
## X <-> X 0.3540660 0.27419454 0.3717394
## X <-> Y 0.7424551 0.03558023 0.2219647
## Y <-> Y 0.7818313 0.06264417 0.1555245
##
## $All_correlations
## r_A r_C r_E Total
## X <-> X 1 1.0000000 1.0000000 1.0000000
## X <-> Y 1 -0.1923834 0.6541744 0.6582171
## Y <-> Y 1 1.0000000 1.0000000 1.0000000# Formatting the table for html or latex report
# Powered by knitr and kableExtra
library(kableExtra)
M %>%
behead() %>%
kable(digits = 3,
align = 'c') %>%
add_complex_header_above(M)|
A(%)
|
C(%)
|
E(%)
|
|||||||
|---|---|---|---|---|---|---|---|---|---|
| est | lCI | uCI | est | lCI | uCI | est | lCI | uCI | |
| X | 0.354 | 0.139 | 0.582 | 0.274 | 0.092 | 0.486 | 0.372 | 0.289 | 0.472 |
| Y | 0.782 | 0.505 | 0.879 | 0.063 | 0.000 | 0.335 | 0.156 | 0.115 | 0.213 |
# Saving the table for the output
M %>%
register_output(
name = 'Table 1',
caption = 'Table 1: Variance components')
## A(%)------------------------- C(%)----------------------------- E(%)-------------------------
## est lCI uCI est lCI uCI est lCI uCI
## X 0.3540660 0.1388926 0.5823833 0.27419454 9.228158e-02 0.4861962 0.3717394 0.2892442 0.4715775
## Y 0.7818313 0.5048375 0.8786609 0.06264417 6.563101e-31 0.3347188 0.1555245 0.1145364 0.2128753